Chapter 7 Diversity
7.1 Alpha diversity
# Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 0) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(richness = 1) %>%
rownames_to_column(var = "sample")
neutral <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(neutral = 1) %>%
rownames_to_column(var = "sample")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, tree = genome_tree) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(phylogenetic = 1) %>%
rownames_to_column(var = "sample")
# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
to.elements(., GIFT_db) %>%
traits2dist(., method = "gower")
functional <- genome_counts_filt %>%
filter(genome %in% labels(dist)[[1]]) %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, dist = dist) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(functional = 1) %>%
rownames_to_column(var = "sample") %>%
mutate(functional = if_else(is.nan(functional), 1, functional))
# Merge all metrics
alpha_div <- richness %>%
full_join(neutral, by = join_by(sample == sample)) %>%
full_join(phylogenetic, by = join_by(sample == sample)) %>%
full_join(functional, by = join_by(sample == sample))alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == sample)) %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Species",
breaks=c("Pk","Eb","Ha"),
labels=c("Pipistrellus kuhlii","Eptesicus bottaew","Hipsugo ariel"),
values=c("#e5bd5b", "#6b7398","#e2815a")) +
scale_fill_manual(name="Species",
breaks=c("Pk","Eb","Ha"),
labels=c("Pipistrellus kuhlii","Eptesicus bottaew","Hipsugo ariel"),
values=c("#e5bd5b50", "#6b739850","#e2815a50")) +
facet_wrap(. ~ metric, scales = "free", ncol=4) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_blank())7.2 Beta diversity
beta_q0n <- genome_counts %>%
select(where(~!all(. == 0))) %>% # remove empty samples
column_to_rownames(., "genome") %>%
hillpair(., q = 0)
beta_q1n <- genome_counts %>%
select(where(~!all(. == 0))) %>% # remove empty samples
column_to_rownames(., "genome") %>%
hillpair(., q = 1)
beta_q1p <- genome_counts %>%
select(where(~!all(. == 0))) %>% # remove empty samples
column_to_rownames(., "genome") %>%
hillpair(., q = 1, tree = genome_tree)
beta_q1f <- genome_counts %>%
select(where(~!all(. == 0))) %>% # remove empty samples
filter(genome %in% labels(dist)[[1]]) %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1, dist = dist)7.2.1 Richness (q0n)
| Df | Sum Sq | Mean Sq | F | N.Perm | Pr(>F) |
|---|---|---|---|---|---|
| 2 | 0.0005478071 | 0.0002739036 | 0.3643863 | 999 | 0.695 |
| 69 | 0.0518662354 | 0.0007516846 | NA | NA | NA |
| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| species | 2 | 0.004372878 | 0.03552136 | 1.270621 | 0.239 |
| Residual | 69 | 0.118732719 | 0.96447864 | NA | NA |
| Total | 71 | 0.123105597 | 1.00000000 | NA | NA |
beta_q0n$C %>%
vegan::metaMDS(., trymax = 500, k = 2, trace = 0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(species) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = species)) +
scale_color_manual(values = c("#e5bd5b", "#6b7398","#e2815a")) +
scale_shape_manual(values = 1:10) +
geom_point(size = 4) +
# stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
)7.2.2 Neutral (q1n)
| Df | Sum Sq | Mean Sq | F | N.Perm | Pr(>F) |
|---|---|---|---|---|---|
| 2 | 0.2034005 | 0.10170027 | 7.084138 | 999 | 0.003 |
| 69 | 0.9905677 | 0.01435605 | NA | NA | NA |
| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| species | 2 | 1.567475 | 0.06477038 | 2.389336 | 0.003 |
| Residual | 69 | 22.633010 | 0.93522962 | NA | NA |
| Total | 71 | 24.200485 | 1.00000000 | NA | NA |
beta_q1n$C %>%
vegan::metaMDS(., trymax = 500, k = 2, trace = 0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(species) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = species)) +
scale_color_manual(values = c("#e5bd5b", "#6b7398","#e2815a")) +
scale_shape_manual(values = 1:10) +
geom_point(size = 4) +
# stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
)7.2.3 Phylogenetic (q1p)
| Df | Sum Sq | Mean Sq | F | N.Perm | Pr(>F) |
|---|---|---|---|---|---|
| 2 | 0.135256 | 0.06762802 | 3.065569 | 999 | 0.057 |
| 69 | 1.522175 | 0.02206051 | NA | NA | NA |
| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| species | 2 | 0.4397545 | 0.04243533 | 1.528898 | 0.139 |
| Residual | 69 | 9.9231780 | 0.95756467 | NA | NA |
| Total | 71 | 10.3629324 | 1.00000000 | NA | NA |
beta_q1p$C %>%
vegan::metaMDS(., trymax = 500, k = 2, trace = 0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(species) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = species)) +
scale_color_manual(values = c("#e5bd5b", "#6b7398","#e2815a")) +
scale_shape_manual(values = 1:10) +
geom_point(size = 4) +
# stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
)7.2.4 Functional (q1f)
| Df | Sum Sq | Mean Sq | F | N.Perm | Pr(>F) |
|---|---|---|---|---|---|
| 2 | 0.2176262 | 0.10881312 | 2.510851 | 999 | 0.098 |
| 69 | 2.9902627 | 0.04333714 | NA | NA | NA |
| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| species | 2 | 0.9379128 | 0.0763021 | 2.849874 | 0.044 |
| Residual | 69 | 11.3541841 | 0.9236979 | NA | NA |
| Total | 71 | 12.2920969 | 1.0000000 | NA | NA |
beta_q1f$C %>%
vegan::metaMDS(., trymax = 500, k = 2, trace = 0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(species) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = species)) +
scale_color_manual(values = c("#e5bd5b", "#6b7398","#e2815a")) +
scale_shape_manual(values = 1:10) +
geom_point(size = 4) +
# stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
)